Orbital Stability of Exomoons and Submoons with Applications to Kepler 1625b-I

In stellar systems, planets orbit stars and moons orbit those planets in an hierarchical configuration (i.e., a satellite chain). Naturally, one questions whether exomoons, or satellites, could host their own satellites (i.e., submoons) and thereby extend the chain one link further. In the solar system, only two planets (Mercury and Venus) are without moons, where there are ∼350 minor planets or asteroids that also have natural satellites.⁴ For the host bodies to hold onto their moons for long timescales, the satellites must reside within a restricted dynamical space (Domingos et al. 2006, hereafter DWY06) that depends on both the Roche and Hill radius as natural boundaries. The ideal conditions for moons occur when the tidal stresses of the surrounding larger objects are small as to not affect the orbital stability of the satellite. Tidal stress can make or break a satellite by altering its spin and orbital parameters (Reid 1973), as well as potentially disrupting the satellite completely. Tidal evolution in the context of star–planet–moon systems has been studied in the general case of moons orbiting exoplanets (Barnes & O'Brien 2002; Sasaki et al. 2012; Sasaki & Barnes 2014; Piro 2018). Depending on the configuration, moons may migrate inward toward the Roche radius to breakup or collide with their host planet, migrate outward toward the Hill radius leading to expulsion from the system, or they might migrate only temporarily (Barnes & O'Brien 2002; Sasaki et al. 2012; Piro 2018). The tidal migration of the moons is also tied to the host planet's rotation, which in turn is affected by stellar tidal friction (Burns 1973; Ward & Reid 1973).

Submoons (i.e., satellites of satellites) appear to be absent within the solar system, but may exist elsewhere in the Milky Way. To detect such objects, a good starting place will be in the search for exoplanets with moons (i.e., exomoons) and leverage the strengths of the current proposed techniques for discovering exomoons. Such methods include microlensing (Han & Han 2002; Han 2008; Liebig & Wambsganss 2010), direct imaging (Cabrera & Schneider 2007; Agol et al. 2015), cyclotron radio emission (Noyola et al. 2014, 2016), photometric transit timing (Sartoretti & Schneider 1999; Kipping 2009a, 2009b), and radial velocity (RV; Vanderburg et al. 2018). The most promising of these methods are through photometric transit timing and cyclotron radio emission. However, either a degeneracy exists within these methods that makes the detection ambiguous or the current state of technology has yet to achieve the necessary precision to obtain an adequate observation.

Currently, there is only one exomoon candidate, Kepler 1625b-I (Teachey & Kipping 2018), whose host planet is Jupiter-like belonging to a system with a Sun-like host star. The initial vetting of this candidate used data obtained from the Kepler mission as well as follow-up observations with the wide-field camera of the Hubble Space Telescope. Teachey & Kipping (2018) examined several different models after the data was corrected for known systematics of the respective instruments aboard the telescopes and concluded that a Neptune-sized satellite represented the best-fit model quite well. Other observers have contested these results suggesting that the exomoon signal is an artifact of the analysis (Kreidberg et al. 2019) or possibly a misidentified event from other exoplanets orbiting the same host star (Heller 2018). Upon inspection of all the data, Heller et al. (2019) found strong statistical evidence using the Bayesian information criterion which favored the planet–moon model over the single-planet model, thus suggesting the existence of an exomoon.

While Hamers & Portegies Zwart (2018) identified tidal capture as an alternate formation pathway for which the Jovian planet could have acquired such a large moon instead forming the satellite from a circumplanetary disk (Canup & Ward 2006), through numerical simulations, the wide orbit of Kepler 1625b-I could be reproduced if the exomoon candidate begins as a circumstellar co-orbital body with the original core of the giant planet Kepler-1625b. The smaller co-orbital body is later pulled down and captured into a circumplanetary orbit, becoming an exomoon (Hansen 2019). Even the habitability of a putative Earth-like exomoon was studied (Williams et al. 1997; Forgan 2018, 2019) given that the Jupiter-like host exists within the habitable zone of its host star. Evidently, additional observations are needed to confirm the exomoon signal but this may also prove difficult depending on the observed geometry of the exomoon's orbit relative to the host star's radius (Martin et al. 2019).

There are challenges for observing exomoons and likely additional hurdles for submoons, but the current exoplanet candidate Kepler 1625b-I provides the means to address potential outcomes numerically. DWY06 produced a fitting formula to identify the outer boundary for stability concerning exomoons that was informed by earlier works examining the stability of similar architectures (Holman & Wiegert 1999). In the past decade, the potential stability of exomoons was probed for individual systems (Quarles et al. 2012; Cuntz et al. 2013) and more general studies (Quarles et al. 2020) similar to those by Holman & Wiegert (1999). The stability of artificial submoons (e.g., satellites orbiting Earth's moon) is also dependent on the mass distribution within the host satellite where mass concentrations, or mascons, can act to perturb the orbit of the submoon and limit its potential stability, especially for low-altitude orbits. The Lunar Prospector Mission hinted at these effects (Konopliv et al. 2001) and the Gravity Recovery and Interior Laboratory (GRAIL) mission showed prominent mascons, the largest of which are associated with nearside basins, as well as the broad structure of the highlands. It is likely that orbits of submoons orbiting terrestrial satellites can be unstable due to perturbations of the host planet as well as perturbations created by mascons in the host satellite. For the purpose of our study, we do not consider the effect of mascons, and thus we assume host planets and moons exert a uniform gravitational field.

In this work, we address whether this formalism can be improved for exomoons or expanded to include submoons. Besides the outer stability limit, Kollmeier & Raymond (2019) showed that not just any submoon would be orbitally stable due to tidal interactions, thus we investigate the long-term stability of submoons. Future space and ground-based telescopes could overcome some of the challenges in the observing for the photometric, RV, and/or cyclotron radio emission methods, for which we identify the thresholds that such methods should meet in order to confirm the possible existence of a submoon.

In Section 2 we outline our numerical methods that include N-body simulations as well as the possible effects on observations. We discuss our results for the stability of exomoons in Section 3.1.1 and submoons in Section 3.1.2. The possible limits due to tidal migration for an exomoon and submoon orbiting Kepler 1625b are explored in Section 3.2. Section 3.3 explores the potential observational scenarios for Kepler 1625 assuming that a stable exomoon and submoon exist there. Finally, we conclude this work with a summary along with the prospects for future observations in Section 4.

We examine the stability of submoons through N-body simulations for both short (10⁵ yr) and long (10 Myr) timescales. In order to investigate the short timescale, we use the IAS15 integrator in REBOUND (Rein & Liu 2012; Rein & Spiegel 2015), which uses an adaptive time step to ensure adequate numerical accuracy within our simulations. For evaluating the orbital evolution on long timescales, we use a modified version of the N-body integrator mercury6 (Chambers et al. 2002) which has been designed to efficiently evolve a similar hierarchy. In mercury6, we adjust the parameters for the switch-over function so that the evolution of the submoon is evaluated via the Bulirsch–Stoer algorithm and the orbits of the other bodies (planet and moon) are evolved through the Wisdom–Holman (WH) portion of the hybrid algorithm. Since stable submoons must orbit their parent bodies closely, this algorithm is well suited for this particular hierarchy. From this division, we choose an initial time step of 5% of the exomoon orbital period such that the WH integration also maintains a high accuracy over long timescales.

This setup also assumes that some formation pathway exists for both the exomoon and submoon (e.g., tidal capture; Hamers & Portegies Zwart 2018), where the submoon can persist past the early stages of formation of the host planet, even though such events are probably rare. Sasaki et al. (2012) showed that a Jupiter-mass planet can host an Earth-mass satellite and such a satellite will be safe from tidal effects on 10 Gyr timescales at planet separations larger than 0.3 au. We follow the results from Kollmeier & Raymond (2019) and prescribe an exomoon radius such that a long-lived submoon will be negligibly affected by the tidal interactions with the host satellite using

$$\begin{eqnarray}&&{R}_{\mathrm{sat}}\geqslant {\left[\displaystyle \frac{39{M}_{\mathrm{sub}}{k}_{2,\mathrm{sat}}T\sqrt{G}}{2{\left(4\pi {\rho }_{\mathrm{sat}}/3\right)}^{8/3}{Q}_{\mathrm{sat}}}{\left(\displaystyle \frac{3{M}_{{\rm{p}}}}{{\left({{fa}}_{\mathrm{sat}}\right)}^{3}}\right)}^{13/6}\right]}^{1/3},\end{eqnarray} \tag{ 1 }$$

which includes the exomoon radius R_sat, submoon mass M_sub, planet mass M_p, exomoon Love number ${k}_{2,\mathrm{sat}}$, exomoon density ${\rho }_{\mathrm{sat}}$, exomoon tidal quality factor Q_sat, exomoon semimajor axis a_sat, and system lifetime T. Following the system parameters from Teachey & Kipping (2018), a Jovian (4 M_J) planet orbits a Sun-like star (1.04 M_⊙) on a ∼287 day (∼0.86 au) orbit. At such a large separation, the impact of tides from the host star is small on the satellites (and submoons) of Kepler 1625b. For the exomoon candidate, Kepler 1625b-I, we consider a 22 day (or a_sat = 0.024 au) orbit around the Jovian host and a Neptune-like mass (m_sat = 0.0564 M_J) for the exomoon, consistent with the observations from Teachey & Kipping (2018). From the calculated radius in Equation (1), most submoons would likely not support a substantial atmosphere so we use an asteroid-like density (3 g cm⁻³) to calculate the mass of the submoon, which is very small (m_sub = 1 × 10⁻⁶ M_⊕) compared to the host planet and exomoon.

Our general simulations for exomoons use our short timescale (10⁵ yr) in a similar manner as DWY06, where a Jupiter-mass planet begins on a Earth-sized orbit (a_p = 1 au) and hosts an Earth-mass exomoon. The orbital eccentricity of the planet and exomoon are varied from 0.0 to 0.5 in steps of 0.01. Both orbits begin coplanar to one another (i_p = i_sat = 0°) with the argument of pericenter and ascending node orbital elements set to zero (ω = Ω = 0°). The planet begins with its mean anomaly equal to zero (MA_p = 0°), while 20 values of the exomoon's mean anomaly are randomly selected from a uniform distribution from 0° to 360°. The simulations that differ in orbital phase are used to calculate the fractional stability f_stab of an initial condition, where a simulation is given either zero weight or full weight depending on the numerical outcome (unstable or stable, respectively). We use the planet's Hill radius (${R}_{{\rm{H}},{\rm{p}}}={a}_{{\rm{p}}}{\left({\mu }_{{\rm{p}}}/3\right)}^{1/3};$ where μ_p = M_p/M_⊙) to scale our exomoon results. Using the prograde results from DWY06, we increment the exomoon semimajor axis a_sat from 0.25 R_H,p to 0.55 R_H,p. A simulation is terminated and marked unstable if the parent–satellite separation r is less than the parent radius (i.e., collision) or greater than the Hill radius of the parent body (i.e., escape).

We perform two sets of simulations for submoons using a Neptune-like exomoon host, where one set is a grid search similar to the method described above and the other explores more general initial orbits (i.e., all angles drawn from distributions). For the first set, we select a_sat for a 22 day orbit (i.e., 0.22 R_H,p) and vary the submoon semimajor axis a_sub from 0.05 R_H,sat to 0.55 R_H,sat using the exomoon's Hill radius (${R}_{{\rm{H}},\mathrm{sat}}={a}_{\mathrm{sat}}{\left({\mu }_{\mathrm{sat}}/3\right)}^{1/3};$ where ${\mu }_{\mathrm{sat}}={M}_{\mathrm{sat}}/{M}_{{\rm{p}}}$). We choose 0.05 R_H,sat as the lower limit for a_sub because it is approximately the Roche limit, assuming a Neptune-like and asteroid-like density for the host exomoon and submoon, respectively. Our second set explores a more limited range for a_sub ranging from 0.3 R_H,sat to 0.4 R_H,sat, uses the system parameters for Kepler 1625b-I, and draws all angles from distributions. Teachey & Kipping (2018) suggested that the exomoon's orbit could be inclined relative to the orbital plane of the planet and thus, we draw the eccentricities and mutual orbital inclinations from Rayleigh distributions, while the rest of the non-fixed orbital parameters are drawn from uniform distributions. Figure 1 illustrates the extent of our initial setup, where Figure 1(c) shows the putative boundaries within the context of the circular restricted three body problem (Eberle et al. 2008; Musielak & Quarles 2014).

Zoom In Zoom Out Reset image size

3.1. Stability Limits of Exomoons and Submoons

An aspect that can be studied without the need of a complete theory of satellite formation is the orbital evolution of satellites since the survivors of such evolution can be potentially observed. In this context, one of the goals of this work entails determining the stability boundary for orbits of possible exomoons and submoons around the exomoon candidate Kepler 1625b-I. The stability boundary is a useful tool because it places a constraint on the maximum orbital period allowed for (sub)satellites, while tidal interactions can limit the minimum orbital period. DWY06 derived expressions of the critical semimajor axis (a_E) in units of the planet's Hill radius R_H,p, where satellites located beyond the critical semimajor axis would escape the gravitational influence of the host planet. More specifically, their work produced an expression for stable prograde and retrograde orbits, where the eccentricities of the planet e_p, and the satellite e_sat are included (see Equations (5) and (6) from Domingos et al. 2006). Retrograde orbits that are typically associated with the irregular satellites of Jupiter (Jewitt & Haghighipour 2007) also have high eccentricities. This is a supremely limiting factor for the stability of exomoons and submoons given the hierarchical nature of our problem and thus retrograde orbits are not within the scope of this work.

3.1.1. Exomoon Stability Limits Revisited

To understand the stability of submoons, first we must address the stability of their host exomoons. Thus, we perform a suite of simulations varying the initial orbital phase MA_sat and the semimajor axis a_sat of the exomoon, as well as the eccentricity e_p of the host planet. We define the quantity f_stab as the fraction of 20 simulations with randomly chosen MA_sat that survive for 10⁵ yr (i.e., Quarles et al. 2018, 2020). Figure 2(a) shows the results of our simulations, where f_stab is color coded and the white cells mark when f_stab = 0. The largest a_sat (with f_stab = 1) for a given e_p marks the stability boundary and is followed by a transition region (colored cells). The black dashed curve marks the best-fit solution from DWY06 for prograde exomoons and the red solid curve represents our best-fit solution to the stability boundary using the Levenberg–Marquardt algorithm in the curve fit function within the scipy module. The light gray curves mark the uncertainty in our estimation. The first two rows of Table 1 provide quantitatively the coefficients for each curve using ${{\boldsymbol{a}}}_{{\bf{crit}}}({R}_{{\rm{H}},{\rm{p}}})={c}_{1}(1-{c}_{2}{e}_{{\rm{p}}})$ as the model. Our results agree with DWY06 when considering a single choice for the orbital phase of the exomoon, but widely disagree when we approach the stability more probabilistically using f_stab. The first coefficient ${c}_{1}\approx 0.4$ (i.e., stability limit for circular orbits) differs by ∼20%, while the second coefficient c₂ for the e_p term differs by only ∼10%. These differences with DWY06 are beyond the uncertainties given and decreases their estimates for the maximum exomoon mass to only 30% of the reported values (DWY06; see their Table 1). Payne et al. (2013) also found c₁ ≈ 0.4 when investigating exomoons within dynamically packed planetary systems, but their results are ambiguous between choaticity and stability through the use of the Lyapnuov exponent.

Zoom In Zoom Out Reset image size

Table 1.  Coefficients for the Critical Semimajor Axis

	${c}_{1}\pm {\sigma }_{1}$	${c}_{2}\pm {\sigma }_{2}$	${c}_{3}\pm {\sigma }_{3}$
Exomoon (DWY06 a)	0.4895 ± 0.0363	1.0305 ± 0.0612	0.2738 ± 0.0240
Exomoon	0.4061 ± 0.0028	1.1257 ± 0.0273	–
Exomoon	0.4031 ± 0.0007	1.1230 ± 0.0041	0.1862 ± 0.0050
Submoon	0.3352 ± 0.0036	0.6642 ± 0.0360	–
Submoon	0.3210 ± 0.0008	0.2757 ± 0.0060	1.0687 ± 0.0050

Notes.

^aThe coefficients ( ${c}_{1}-{c}_{3}$) and the uncertainties ( ${\sigma }_{1}-{\sigma }_{3}$) from Domingos et al. (2006) are listed using the fitting formula for the external boundary, ${{\boldsymbol{a}}}_{{\bf{crit}}}({R}_{{\rm{H}}})={c}_{1}(1-{c}_{2}{e}_{p}-{c}_{3}{e}_{\mathrm{sat}}$). The remaining rows contain coefficients and uncertainties from our work for the same fitting formulas. The dash symbol (–) denotes when a coefficient is not applicable because the respective eccentricity is fixed at zero.

Download table as:  ASCIITypeset image

Figure 2(b) demonstrates our results when the exomoon's eccentricity e_sat also varies and allows us to define a new model, ${{\boldsymbol{a}}}_{{\bf{crit}}}({R}_{{\rm{H}},{\rm{p}}})={c}_{1}(1-{c}_{2}{e}_{{\rm{p}}}-{c}_{3}{e}_{\mathrm{sat}}),$ similar to DWY06. The color coded cells mark the largest a_sat (with f_stab = 1) for a given set of initial eccentricities e_sat and e_p. As one might expect, the dependence on e_p is stronger than for e_sat. This is shown quantitatively in third row of Table 1, where c₂/c₃ ≈ 6. Additionally, we find that the second and third rows agree (within the uncertainties), which provides a sanity check for our numerical results. We provide a fitting formula that better matches more robust numerical simulations for exomoons, but such a formula is reliable in characterizing a population of exomoons and can be inaccurate when describing particular systems. As a result, we suggest a different approach using a lookup table or interpolation map that relies less on statistical averaging, which occurs in deriving a single fitting formula.

3.1.2. Stability Limits for Submoons

Kollmeier & Raymond (2019) previously studied the possible tidal evolution and stability of submoons following the analytical treatment for the solar system moons (Murray & Dermott 1999) and exomoons (Barnes & O'Brien 2002; Sasaki et al. 2012). However, part of their analysis depends on the stability limit for submoons, where small changes in f (Equation (1)) are amplified by the large exponent. We acknowledge that high uncertainties remain in other terms (e.g., ${k}_{2,\mathrm{sat}}$ or Q_sat), but changes in those values result in a much smaller change to R_sat due to the smaller exponent.

To determine the stability limit, we follow a similar methodology as our investigation of exomoon stability (see Section 3.1.1). The major differences are that the submoon semimajor axis is incremented in units of the host satellite's Hill radius R_H,sat, the submoon's initial orbital phase MA_sub is chosen randomly while the host satellite's mean anomaly is not (${\mathrm{MA}}_{\mathrm{sat}}=0^\circ$), the host satellite's mass is set to nearly equal to Neptune (${M}_{\mathrm{sat}}\approx 18{M}_{\oplus }$), and the host satellite's semimajor axis is no longer varied (${a}_{\mathrm{sat}}\approx 0.015$ au). Figure 2(c) demonstrates the stability of a submoon when only the planetary eccentricity and submoon semimajor axis are varied (identical color code as in Figure 2(a)). Highly stable initial conditions (black cells) cluster first near ${a}_{\mathrm{sub}}\sim {R}_{{\rm{H}},\mathrm{sat}}/3$ for low planetary eccentricity (${e}_{{\rm{p}}}\lesssim 0.2$) and then more widely at ${a}_{\mathrm{sub}}\sim {R}_{{\rm{H}},\mathrm{sat}}/4$. Due to the hierarchical nature of these systems, overlap of mean motion resonances can destabilize submoons in a manner similar to planets in binary systems (Mudryk & Wu 2006; Quarles et al. 2020). Hence a gap forms when the period ratio of the satellite to the submoon equals an integer (${P}_{\mathrm{sat}}/{P}_{\mathrm{sub}}=\sqrt{3/{f}^{3}}$) and f is a fraction of the host satellite's Hill radius. The fourth row of Table 1 shows the results of a linear fit, where the first coefficient ${c}_{1}={R}_{{\rm{H}},\mathrm{sat}}/3$ is the most meaningful.

Figure 2(d) demonstrates how the stability limit varies with changes in both the planetary and satellite eccentricity. The orange region for low eccentricity corresponds to conditions where the stability island near ${a}_{\mathrm{sub}}\sim {R}_{{\rm{H}},\mathrm{sat}}/3$ can persist, where most choices of eccentricity push the stability limit to lower values (green to blue region). Eccentricity values near 0.5 (gray region in Figure 2(d)) are wholly unstable because the stability limit lies within the Roche limit of the host satellite. The final row of Table 1 provides coefficients for a fitting function, but such a model employs some averaging. Specifically, the fitting function removes the plateau at low eccentricity (orange region) and smooths out the transition from 0.25 R_H,sat to 0.20 R_H,sat at high eccentricity of the satellite (light green to blue region) in Figure 2. Comparing the fourth and last row in Table 1, there is a significant difference and the same sanity check that we employ for exomoons does not hold. The difference arises because the host satellite's eccentricity is now allowed to vary and it plays a significant role in the stability of its submoons. By default, we would expect that c₂ in the exomoon case should map onto c₃ for submoon stability. This is mostly true, where some dilution occurs due to the main plateau at low eccentricities and the host planet still has some influence on the submoon.

Kollmeier & Raymond (2019) assumed that the stability limit for submoons would largely match that of exomoons when scaled by the appropriate Hill radius. This is clearly not the case as their results underestimate the true limits in the host satellite's radius. Additionally, our Equation (1) differs from Kollmeier & Raymond by a factor of one-third in the denominator that was due to a typographical error in the manuscript (S. Raymond 2020, private communication). If we use our stability limit for f (assuming nearly circular orbits), the minimum satellite radii reported are increased by a factor k $(={(0.4895/0.33)}^{13/2}\sim 13$). Relaxing the assumption for nearly circular orbits to allow for eccentric orbits for the host planet and satellite, the denominator of the correction factor k can be smaller, which makes the importance of the stability limit even greater.

3.1.3. Long-term Stability of Submoons

Our study of the stability of submoons is inherently limited by computational resources, where some may want to see billion year simulations to decide whether an initial condition is stable for the long term. However, a different approach is to use an analytic formula to ensure stability against tidal migration (i.e., Barnes & O'Brien 2002) and perform N-body simulations for roughly a billion orbits of the submoon. Submoons that survive this timescale are likely to be extant in the face of either tides or gravitational perturbations.

We focus on the region between 0.3 R_H,sat to 0.4 R_H,sat of the host satellite and perform simulations for 10 Myr, or two billion orbits of the submoon. The remaining initial conditions for the submoon are chosen randomly, where e_sub ≲ 0.2 and the inclination relative to the host satellite's orbital plane is ≲15° (see Section 2). The initial semimajor axis for the host planet and satellite follow the parameters from Kepler 1625b, but the remaining orbital elements are chosen randomly. We evolved ∼1500 initial conditions for this longer simulation timescale and find that the stable submoons cluster around ${a}_{\mathrm{sub}}\approx 0.33{R}_{{\rm{H}},\mathrm{sat}}$ with low eccentricity ${e}_{\mathrm{sub}}\lesssim 0.05$. The inclination differences between the orbital planes of the planet, satellite, and submoon have a negligible effect on the long-term stability when they are significantly below 40°.

3.2. Parameter Limits for Exomoons and Submoons in Kepler 1625

The tidal migration of moons depends on the assumed value for the tidal Love number k₂ and the tidal quality factor Q of the host body. The tidal Love numbers for solar system bodies have largely been estimated (Lainey 2016) so that relatively accurate analogs can be used for our purposes. However, the tidal quality factor Q is very uncertain and thus, we are forced to vary this parameter over a huge range (10²–10⁵). Despite this uncertainty, we can identify constraints for the exomoon candidate Kepler 1625b-I and explore the limiting mass of a putative submoon.

Figure 3 shows the upper limits for the (a) exomoon mass, (b) host planet radius, (c) submoon mass, and (d) exomoon radius up to respective stability limit for coplanar circular orbits over a 9 Gyr timescale. Each set of curves uses the system parameters from Teachey & Kipping (2018) to define the properties of the host body and also indicates the change in the upper limit with an assumed value for Q (color coded). Figures 3(a) and (b) have a star symbol that denotes the location of Kepler 1625b and its candidate moon within the parameter space. If we assume that Kepler 1625b is Jupiter-like in its tidal quality, then values below (or to the right of) the red dashed curve in Figures 3(a) and (b) are allowed. The exomoon parameters from Teachey & Kipping (2018) suggest that the candidate exomoon is similar to Neptune and thus values below the green dashed curves in Figures 3(c) and (d) are allowed. The maximum submoon mass is 3 × 10⁻⁵ M_⊕, which is slightly less massive than the asteroid Vesta. The minimum exomoon and submoon separations are 0.14 R_H,p and 0.2 R_H,sat, respectively, given the measured radius of each host body and reasonable value of Q (red or green dashed curve, respectively).

Zoom In Zoom Out Reset image size

3.3. Detection Thresholds for Submoons

Although Kollmeier & Raymond (2019) introduced a formalism for the stability of a planet–moon–submoon, the possibility of observing such subsatellites has not been addressed. The observation and confirmation of exomoons has proved challenging through current proposed techniques, as mentioned in Section 1. The exploration of exomoon and submoon detection is necessary to contextualize our work within existing observational capabilities. In this section, we use current exomoon detection techniques and determine the noise thresholds necessary for a detection of a submoon. In particular, we use the transit and RV methods, which have been prolific in exoplanet confirmations. Additionally we discuss the challenges that submoon detection in radio wavelengths poses.

3.3.1. Transit Method

To explore the feasibility of detecting submoons through transits we utilize Python packages REBOUND and batman. We use a sample from our long-term (10 Myr) integrations to produce the synthetic light curves (see Table 2). The package batman (Kreidberg 2015) computes transit light curves for a single exoplanet using analytic expressions introduced by Mandel & Agol (2002) and include quadratic or nonlinear limb darkening. We use the relative locations of each body on the sky plane in our rebound simulations to determine how much area covers the stellar disk and approximate the net flux (exoplanet+exomoon) using a single planet with an equivalent area (i.e., homeomorphic area). We account for partial eclipses, or blends, by finding the intersectional area between two disks on the sky (e.g., exoplanet+star or exoplanet+exomoon) to correct the equivalent area for possible double counting.

Table 2.  Initial Orbital Elements for the Synthetic Light Curve and RV of Kepler 1625

	M (${M}_{\oplus }$)	R (${R}_{\oplus }$)	a (au)	RH (au)	e	i (deg)	ω (deg)	ν (deg)
Exoplanet	1270	10	0.86348384	0.09233077	0.01147105	0	250.31211	190.91273
Exomoon	17.91	4	0.02401178	0.00402236	0.01159945	25.805015	264.07880	322.26965
Submoon	10−6	0.0196	0.00144350	0.00000383	0.09984839	0.92405241	189.40435	22.038811

Note. The semimajor axis a, Hill Radius R_H, eccentricity e, inclination i, argument of periastron ω, and true anomaly ν are given relative to the orbital plane of the respective host body, where we set the ascending node ${\rm{\Omega }}=180^\circ$ for all bodies. After a coordinate transformation to a common frame, we rotate the system by 90° so that the line of sight is along the z-axis.

Download table as:  ASCIITypeset image

We simulate the light curves over a four year period (similar to the Kepler mission), identify four planetary transits around the host star, and evaluate the contribution of the satellites (exomoon and submoon) to the transit depth. The synthetic light curves in Figure 4 show satellite transits in addition to planetary transits, where Figures 4(a) and (d) illustrate when the satellites transit ±1 day relative to the planetary transit t_o.

Zoom In Zoom Out Reset image size

The illustrations in Figure 4 are annotated (I, II, or III) to identify the relative location of the satellites. Due to stability constraints on the submoon (see Section 3.1.2), the exomoon–submoon separation is always less than a stellar radius (${a}_{\mathrm{sub}}\approx {R}_{\star }/8;$ ${R}_{\star }=1.73{R}_{\odot }$) which implies that submoons, in general, will always transit as long as their host satellite transits too. The 125 km submoon produces a transit depth of 0.01 ppm and is much smaller than the expected photometric noise in 30 minute cadence observations. Figure 4(b) shows only a planetary transit because the exomoon separation and relative inclination are large enough for transits to not occur at every planetary transit. Figure 4(c) illustrates the difference in velocity on the sky plane between the planet and satellites, where the satellites transit, followed by the planet+satellites, and ending only with the planet.

3.3.2. RV Method

We study the effect of satellites (exomoon and submoon) on potential RV measurements of Kepler 1625 using REBOUND to produce synthetic radial velocities that consider the reflex motion of either the star or the planet. Traditional RV measurements quantify the reflex motion of a star due to a planet (see Lovis & Fischer 2010 for a review), but Vanderburg et al. (2018) recently showed that RV measurements taken in a direct imaging campaign could, in principle, detect the reflex motion of an exoplanet due to an exomoon. We perform four separate numerical simulations over a single planetary period: (1) star–planet, (2) planet–satellite, (3) satellite–submoon, and (4) planet–satellite–submoon so that we can identify the reflex motion using the line-of-sight velocity v_z of a body directly.

The stellar reflex motion due to a planet ${v}_{\star ,{\rm{p}}}$ in Figure 5(a) shows an RV semi-amplitude of ∼120 m s⁻¹ (red curve), where the inclusion of a large satellite ${v}_{\star ,\{{\rm{p}},\mathrm{sat}}\}$ (black) produces similar results. The exomoon produces a slight shift (∼8 m s⁻¹) in the RV curve indirectly by causing the planet to accelerate coherently with its motion induced by the host star. Assuming optimistic observational conditions, this difference in RV is detectable using current-generation spectrographs. The impact of a submoon on the stellar reflex motion is negligible because the submoon causes a ∼0.3 mm s⁻¹ oscillation in the host satellite's velocity.

Zoom In Zoom Out Reset image size

As Vanderburg et al. (2018) showed, direct imaging observations could uncover exomoons by obtaining RV measurements of the host exoplanet. Figure 5(b) demonstrates the expected RV curve from such measurements, where the planetary reflex velocity ${v}_{{\rm{p}},\mathrm{sat}}$ (blue) has a semi-amplitude (∼150 m s⁻¹) slightly larger the stellar RV from the host planet, ${v}_{\star ,{\rm{p}}}$. This is possible because the mass ratio for the planet–satellite is larger than for the star–planet and the exomoon is much closer to its host body. However, RV measurements of the exoplanet directly will also include the center of mass velocity v_p added to the planetary reflex velocity (i.e., ${v}_{{\rm{p}}}+{v}_{{\rm{p}},\mathrm{sat}}$) and is represented through the black dashed curve in Figure 5(b). Isolating the exomoon signal would require ∼0.5% precision in the RV measurements, which Vanderburg et al. (2018) show is possible with an 8 m class telescope, high-resolution spectrograph, and modern adaptive optics systems. Our submoon produces a negligible effect on the planetary RV, where even the most massive submoons allowed by tidal migration would still produce minuscule shifts in the RV curve and is too small to be identified with current and future instruments.

Although our work arrives at the same general conclusion as Vanderburg et al. (2018), we do note that the most extreme RV semi-amplitudes presented by Vanderburg et al. (2018) are not likely realized because tidal migration (see Figure 3(a)) over the lifetime of the system would have caused the exomoon to: (1) collapse onto the host planet, (2) migrate outward past the stability boundary, or (3) migrate outward to produce smaller RV semi-amplitudes.

3.3.3. Radio Emission

Another detection method takes advantage of the radio signal that an exomoon like the candidate Kepler 1625b-I would emit due to the interaction of a satellite with its host planet's magnetic field. Following Noyola et al. (2014), we calculate intrinsic power P_s and incident flux S of such a radio emission for the exomoon candidate orbiting Kepler 1625b.

The expected power for Kepler 1625b-I will depend on several parameters, such as the magnetic field strength B_s, the plasma speed V₀, and the plasma density ρ_s, where each depends on the separation of the exomoon from its host planet. Noyola et al. (2014) showed the detectability of Io-like exomoons orbiting a Jupiter analog, where we can augment their calculation (Noyola et al. 2014, see their Equation (2)) through appropriate scale factors. Io tightly orbits Jupiter at only 0.008 R_H,p, where the intrinsic power P_s is 5 GW. Signals from such a system with S ≈ 50 μJy are detectable at 50 MHz and only 15 light years away. Kepler 1625b-I is more widely separated (0.26 R_H,p) from its host planet and is much larger (4 R_⊕), which yields a power of 3.2 GW assuming an Io-like plasma density, which is slightly less than Io. The power would rise to 150 GW, if the exomoon orbits at the inner limit (0.12 R_H,p) set by tidal migration given a Jupiter-like dissipation (see Figure 3(a)). However the incident flux is also inversely proportional to the distance d squared to the system (i.e., $S\propto {P}_{s}/{d}^{2}$) and Kepler 1625 is approximately 8000 light years away. As a result, the highest expected flux from the Kepler 1625 system is 10⁻⁴ times smaller than a Jupiter-Io analog, or ∼5 nJy. Applying our values in Table 2, the estimated flux drops even lower to 10 pJy due to a much weaker magnetic field at 0.26 R_H,p.

We find that regardless of the intrinsic power P_s of the planet–moon interaction, its distance from Earth poses the biggest challenge to obtain a detection. We refrain from exploring the effect of a submoon in the radio emission of the moon, as its radius is too small (125 km) to contribute to the intrinsic power. Even if the submooon could boost the power by an order of magnitude through plasma sharing between satellites (Noyola et al. 2016), the incident flux would still remain incredibly small due to the large distance to the system.

The orbital stability of exomoons (satellites) and submoons (subsatellites) largely depends on the orbital eccentricity of the respective host body and is limited to a fraction of the host body's Hill radius, R_H. Although, a larger dependence on the (sub-)satellite's eccentricity may arise in multi-body systems (e.g., Galilean moons). We perform N-body simulations for 10⁵ planetary orbits and use 20 random initial orbital phases to determine a stability limit (in R_H), or the largest, stable semimajor axis, of exomoon and submoon systems. As a result, the stability limit for an exomoon is 40% of the host planet's Hill radius (R_H,p) and for a submoon is 33% of the host satellite's Hill radius (R_H,sat), assuming circular, coplanar orbits. Additionally, we determine constraints on the physical properties of exomoons or submoons so that they are extant over a wide range of tidal migration scenarios. Exomoons are detectable given current photometric capabilities (e.g., Kepler or the Transiting Exoplanet Survey Satellite), where large submoons (∼400 km) require less than 0.3 ppm photometric precision for a Sun-like star. Both exomoons and submoons will require next-generation RV facilities, where the RV amplitude on the host star for exomoons is small and for submoons is negligible. The flux through radio emission through interactions of Kepler 1625b-I on the planet's magnetic field would be incredibly small, $[{10}^{-9}-{10}^{-12}]$ Jy, and would require extremely long integration times ($\gt {10}^{3}$ hr) with a telescope like the Square Kilometer Array. However, if the system were closer (∼4.6 pc) the emission would be detectable since the incident flux would be much larger and would also require significantly lower integration time.

DWY06 performed the quintessential study of exomoon stability. However they probed only a single trajectory (orbital phase) for a given pair of initial conditions (${a}_{\mathrm{sat}},{e}_{\mathrm{sat}}$) and caused them to overestimate the outer boundary ${{\boldsymbol{a}}}_{{\bf{crit}}}$ (a_E in their notation) of stability by 20%. Our simulations are evolved for 10× more planetary orbits and consider 20 trajectories per initial condition (see Figure 2 ⁵ ). Using these simulations, we update the coefficients for the fitting formula used by DWY06 (see Table 1). In addition, we expand upon the methodology to consider the stability limit for submoons (in R_H,sat) and find that submoons are stable to only ${R}_{{\rm{H}},\mathrm{sat}}/3$ (see Table 1).

Analytical tidal migration estimates from Barnes & O'Brien (2002) depend on prior estimates of the stability limit as a fraction f of the host planet's Hill radius and are typically used to identify the limiting physical characteristics of exomoons (DWY06). These expressions raise f to a very high exponent (6.5), where large deviations can result from small changes. The major difference for exomoons is that large satellite masses with ($0.4{R}_{{\rm{H}},{\rm{p}}}\lt {a}_{\mathrm{sat}}\leqslant 0.49{R}_{{\rm{H}},{\rm{p}}}$) are not viable due to tidal migration. Kollmeier & Raymond (2019) based their analysis for submoons on Barnes & O'Brien (2002) and DWY06, where the differences for submoons can be much larger due to a less smooth stability surface (see Figure 2). We use Equation (8) from Barnes & O'Brien (2002) and the exomoon candidate system Kepler 1625b-I (Teachey & Kipping 2018) to determine the upper mass limit (70% the mass of Vesta) and upper radius limit (375 km) for submoons. We also show that a Neptune-like exomoon would be orbitally stable around Kepler 1625b for most reasonable assumptions on the tidal quality (Q ≳ 10³) of the host planet.

We evaluate the detection prospects of exomoons and submoons orbiting Kepler 1625b through a range of methods, such as those using photometry, spectroscopy, and radio emission. In the photometric method, the transit depth of the exomoon candidate (Kepler 1625b-I) is distinguishable from the host planet transits, but a submoon would require an extreme photometric precision below 1 ppm. We produce synthetic RVs that show the stellar reflex motion or the planetary reflex motion because the latter could be measured using large telescopes with state-of-the-art direct imaging capabilities. Through both of these measures, we find that the presence of exomoons is identifiable, where an 8 m s⁻¹ difference arises in the stellar RV and 150 m s⁻¹ oscillation in the planetary RV. Recent observations have uncovered star–planet interactions through radio emissions (Vedantham et al. 2020), but such techniques require a much higher sensitivity (nJy scale) to make a detection of an exomoon orbiting Kepler 1625b.

It is worth noting that these observational techniques carry their own biases and in addition to stability, these biases may affect the discovery of exomoons and submoons. The transit method is biased toward well-separated large (sub-)satellites, the RV method has a bias toward more massive bodies, and direct imaging includes both of these biases. Our study for exomoons and submoons shows that orbital stability favors these biases and it makes sense that the first exomoon candidate by Teachey & Kipping (2018), Kepler 1625b-I, is large and massive. However, these features might not be common in the general population of exomoons or submoons.

The authors thank Sean Raymond for his thorough review with comments that enhanced the quality of the manuscript. We are grateful to Jason Barnes for useful comments and discussion. M.R.F acknowledges support from the NRAO Gröte Reber Fellowship and the Louis Stokes Alliance for Minority Participation Bridge Program at the University of Texas at Arlington. This research was supported in part through research cyberinfrastructure resources and services provided by the Partnership for an Advanced Computing Environment (PACE) at the Georgia Institute of Technology.

Software: mercury6 (Chambers et al. 2002), batman (Kreidberg 2015), rebound (Rein & Liu 2012; Rein & Spiegel 2015).